#include "BSpline.hpp"
#include "CubicSpline.hpp"
#include "LinearSpline.hpp"
#include "SplineBase.hpp"
#include "matplotlibcpp.h"
#include <cmath>
#include <iostream>
#include <vector>

using namespace std;
namespace plt = matplotlibcpp;

const double PI = 3.1415926;

void TestLinearSpline() {
  // 实现了线性插值
  vector<double> x = {0, 1, 1.5, 3, 5};
  vector<double> y = {0, 3, -1, 2, 3};
  BSpline bspline(1, x, y);
  LinearSpline linspline(x, y);
  vector<double> xData;
  vector<double> yb;
  vector<double> ypp;
  int num = 1000;
  for (int i = 0; i <= num; i++) {
    xData.push_back(5.0 * i / num);
    yb.push_back(bspline.interpolate(5.0 * i / num));
    ypp.push_back(linspline.interpolate(5.0 * i / num));
  }
  plt::figure_size(1200, 700);
  plt::plot(xData, ypp, "b");
  plt::plot(xData, yb, "r--");
  plt::save("../figure/Linear.png", 300);
}

void CheckPPEqualB() {
  // 测试两种三次样条拟合sinx，验证在三种边界条件下得到的样条是一样的
  int num_of_points = 10;
  vector<double> xData;
  vector<double> yData;
  for (int i = 0; i <= num_of_points; i++) {
    double x = 2 * PI * i / num_of_points;
    xData.push_back(x);
    yData.push_back(sin(x));
  }

  CubicSpline ppform(xData, yData, SplineBase::natural);
  BSpline bform(3, xData, yData, BSpline::natural);

  int num = 1000;
  vector<double> x(num + 1);
  vector<double> ypp(num + 1);
  vector<double> yb(num + 1);
  for (int i = 0; i < num + 1; i++) {
    x[i] = i * 2 * PI / num;
    ypp[i] = ppform.interpolate(x[i]);
    yb[i] = bform.interpolate(x[i]);
  }

  plt::figure_size(1200, 700);
  plt::plot(x, ypp, "b");
  plt::plot(x, yb, "r--");
  plt::save("../figure/natural.png", 300);

  ppform.setBoundary(SplineBase::complete, 1, 1);
  bform.setBoundary(BSpline::complete, 1, 1);
  ppform.setupSplines();
  bform.setupSpline();

  for (int i = 0; i < num + 1; i++) {
    ypp[i] = ppform.interpolate(x[i]);
    yb[i] = bform.interpolate(x[i]);
  }

  plt::figure_size(1200, 700);
  plt::plot(x, ypp, "b");
  plt::plot(x, yb, "r--");
  plt::save("../figure/complete.png", 300);

  ppform.setBoundary(SplineBase::periodic);
  bform.setBoundary(BSpline::periodic);
  ppform.setupSplines();
  bform.setupSpline();

  for (int i = 0; i < num + 1; i++) {
    ypp[i] = ppform.interpolate(x[i]);
    yb[i] = bform.interpolate(x[i]);
  }

  plt::figure_size(1200, 700);
  plt::plot(x, ypp, "b");
  plt::plot(x, yb, "r--");
  plt::save("../figure/periodic.png", 300);
}

void MoreTest() {
  // 测试函数 sin(1/x)
  int num_of_points = 100;
  vector<double> xData;
  vector<double> yData;
  double start = 0.02;
  double end = 0.2;
  for (int i = 0; i <= num_of_points; i++) {
    double x = end * i / num_of_points + start;
    xData.push_back(x);
    yData.push_back(sin(1.0 / x));
  }

  BSpline bform(3, xData, yData);

  int num = 1000;
  vector<double> x(num + 1);
  vector<double> y(num + 1);
  vector<double> yb(num + 1);
  for (int i = 0; i < num + 1; i++) {
    x[i] = end * i / num + start;
    y[i] = sin(1.0 / x[i]);
    yb[i] = bform.interpolate(x[i]);
  }

  plt::figure_size(1200, 700);
  plt::plot(x, yb, "b");
  plt::plot(x, y, "r--");
  plt::save("../figure/test.png", 300);

  // 测试函数 sinx/x
  num_of_points = 19;
  start = -20;
  end = 20;
  xData.clear();
  yData.clear();
  for (int i = 0; i <= num_of_points; i++) {
    double x = (end - start) * i / num_of_points + start;
    xData.push_back(x);
    yData.push_back(10.0 * sin(x) / x);
  }

  BSpline bform_2(3,xData,yData);

  for (int i = 0; i < num + 1; i++) {
    x[i] = (end - start) * i / num + start;
    y[i] = 10.0 * sin(x[i]) / x[i];
    yb[i] = bform_2.interpolate(x[i]);
  }

  plt::figure_size(1200, 700);
  plt::plot(x, yb, "b");
  plt::plot(x, y, "r--");
  plt::save("../figure/test2.png", 300);
}

void AnyDegree(){
  // 测试任意次样条，使用函数 exp(-x^2)
  int num_of_points = 10;
  vector<double> xData;
  vector<double> yData;
  double start = -2;
  double end = 2;
  for (int i = 0; i <= num_of_points; i++) {
    double x = (end - start) * i / num_of_points + start;
    xData.push_back(x);
    yData.push_back(exp(-x * x));
  }

  int num = 1000;
  vector<double> x(num + 1);
  vector<double> y(num + 1);
  vector<double> y4(num + 1);
  vector<double> y5(num + 1);
  BSpline bform_4(7, xData, yData);
  BSpline bform_5(5, xData, yData);
  for (int i=0; i<=num; i++) {
    x[i] = (end - start) * i / num + start;
    y[i] = exp(-x[i] * x[i]);
    y4[i] = bform_4.interpolate(x[i]);
    y5[i] = bform_5.interpolate(x[i]);
  }
  plt::figure_size(1200, 700);
  plt::plot(x, y4, "b");
  plt::plot(x, y5, "r");
  plt::plot(xData, yData, "o");
  plt::save("../figure/anydegree.png", 300);
}

int main() {
  std::cout << "----- running More Test -----" << std::endl;
  TestLinearSpline();
  CheckPPEqualB();
  MoreTest();
  AnyDegree();
  std::cout << "----- complete Test -----" << std::endl;
}